f <- function(x)
  {
    answer <- x*x*x
    return(answer)
  }

Sad_X <- function (x, y, h, uniform) {
  # Dviem atsitiktiniais dydziais generuojame X_h (log(S_h))
  
  x1 = x - sqrt(y*h)
  x2 = x + sqrt(y*h)
  
  p = 1/2
    
  if (uniform < p) {
    answer <- x1
  }
  else {
    answer <- x2
  }
  
  return (answer)
}


Sad_rho_e <- function (x, y, h, rho, sigma, X_h, Y_h) {
  # Dviem atsitiktiniais dydziais generuojame X_h
  # Rho <> 0
  
  #tarpinis = sqrt(1-rho*rho) * (S_h-s) + 
  #  rho * (s/sigma) * (Y_h-y)
  
  #  if (tarpinis < 0) 
  #    print ("Neigiamas Sad_rho tarpinis")
  
  answer = 
    x + 
    sqrt(1-rho*rho) * (X_h-x) + 
    rho * (1/sigma) * (Y_h-y)

  #exp_answer = exp (answer)
  
  return (answer)
}


Sad_rho <- function (s, y, h, rho, sigma, S_h, Y_h) {
# Dviem atsitiktiniais dydziais generuojame S_h
# Rho <> 0

  tarpinis = sqrt(1-rho*rho) * (S_h-s) + 
    rho * (s/sigma) * (Y_h-y)
  
#  if (tarpinis < 0) 
#    print ("Neigiamas Sad_rho tarpinis")
  
	answer = 
		s + 
		sqrt(1-rho*rho) * (S_h-s) + 
		rho * (s/sigma) * (Y_h-y)

	if (answer < 0) 
		print ("Neigiamas Sad_rho")

	return (answer)
}

Sad <- function (s, y, h, uniform) {
# Dviem atsitiktiniais dydziais generuojame S_h

	x1 = s + y*s*h - sqrt((y+y*y*h)*s*s*h)
	x2 = s + y*s*h + sqrt((y+y*y*h)*s*s*h)

	p = s/(2*x1)

	if (x1 < 0)
		print ("Neigiamas x1 Sad")
	if (x2 < 0)
		print ("Neigiamas x2 Sad")

	if (uniform < p) {
		answer <- x1
	}
	else {
		answer <- x2
	}

	return (answer)
}

Sad4 <- function (s, y, h, uniform) {
# Keturiais atsitiktiniais dydziais generuojame S_h

	m1 = s
	m2 = 4*h*s*s*y*(1/4 + y*h) + s*s
	m3 = 1/3*y^2*s^3*h^2 + s^3 + 12*h*s^3*y*(1/4+y*h)
	m4 = 5*y^2*s^4*h^2+7*s^4+56*h*s^4*y*(1/4 + y*h)

	es = (m1*m4-m2*m3)/(m1*m3-m2^2)
	ku = (m2*m4-m3^2)/(m1*m3-m2^2)

	delta = es^2 - 4*ku

	x1 = (es^2-sqrt(delta))/2
	x2 = (es^2+sqrt(delta))/2

	p1 = (m1*x2-m2)/(x1*(x2-x1))
	p2 = (m2-m1*x1)/(x2*(x2-x1))

	print ("Sad4 ")
	print (s)

	if (uniform < p1) {
		answer <- x1
	}
	else 
	if (uniform < (p1+p2)) {
		answer <- x2
	}
	else {
		answer <- 0
	}

	return (answer)
}


Sad4c <- function (s, y, h, uniform) {
# Keturiais atsitiktiniais dydziais generuojame S_h centruotu budu

	m1 = 0
	m2 = 4*h*s^2*y*(1/4 + y*h)
	m3 = 1/3*y^2*s^3*h^2
	m4 = 3*y^2*s^4*h^2

	es = (m1*m4-m2*m3)/(m1*m3-m2^2)
	ku = (m2*m4-m3^2)/(m1*m3-m2^2)

	delta = es^2 - 4*ku

	x1 = (es^2-sqrt(delta))/2
	x2 = (es^2+sqrt(delta))/2

	p1 = (m1*x2-m2)/(x1*(x2-x1))
	p2 = (m2-m1*x1)/(x2*(x2-x1))

	print ("Sad4 ")
	print (s)

	print (x1)
	print (x2)

	print (p1)
	print (p2)


	if (uniform < p1) {
		answer <- (s+x1)
	}
	else 
	if (uniform < (p1+p2)) {
		answer <- (s+x2)
	}
	else {
		answer <- s
	}

	return (answer)
}


CIRad <- function (y0, h, a, uniform) {
# Dviem atsitiktiniais dydziais generuojame Y_h (pagal CIR)
# a = sigma*sigma
	
	x1 = y0 + a*h - sqrt((y0+a*h)*a*h)
	x2 = y0 + a*h + sqrt((y0+a*h)*a*h)

	p = y0/(2*x1)

	if (uniform < p) {
		answer <- x1
	}
	else {
		answer <- x2
	}

	return (answer)
}

CIRad4 <- function (y, h, a, uniform) {
# Keturiais atsitiktiniais dydziais generuojame Y_h (pagal CIR)
# a = sigma*sigma

	m1 = y
	m2 = y^2 + y*a*h

	es = (4*y^2+9*y*a*h + 3*(a*h)^2)/(2*y+a*h)
	delta = (a*h*(16*y^3 + 33*y^2*a*h + 18*y*a^2*h^2 + 3*a^3*h^3))/(2*y + a*h)^2
	
	x1 = (es^2-sqrt(delta))/2
	x2 = (es^2+sqrt(delta))/2

	p1 = (m1*x2-m2)/(x1*(x2-x1))
	p2 = (m2-m1*x1)/(x2*(x2-x1))

	print ("CIRad4 ")
	print (y)
	print (p1)
	print (p2)
	print (uniform)

	if (uniform < p1) {
		answer <- x1
	}
	else 
	if (uniform < (p1+p2)) {
		answer <- x2
	}
	else {
		answer <- 0
	}

	return (answer)
}


ad <- function (s, y0, h, sigma, uniform) {
# Pirmas bandymas dviem atsitiktiniais dydziais generuoti S_h...

	#if (y > 0) {

		m1 = s
		m2 = s*s* exp ( ( (sqrt(2)*y0)/sigma ) * tan ( (sigma*h)/sqrt(2) ) )
		#print (m2)

		x1 = (m2 - sqrt (m2*m2-m1*m1*m2)) / m1

		x2 = (m2 + sqrt (m2*m2-m1*m1*m2)) / m1

		#print (y)
		#print (x1)
		#print (x2)
		p = m1 / (2*x1)
		#print (uniform)
		#print (p)

		if (uniform < p) {
			 answer <- x1
			#print ("x1")
		}
		else {
			answer <- x2
			#print ("x2")
		}
	#}
	#else {
	#	answer = 0
	#}
	
	return (answer)
}

ad2 <- function (s, y0, h, sigma, uniform) {
# Senas bandymas 4 atsitiktiniais dydziais generuoti S_h

	m1 = s
	m2 = s*s*exp(((sqrt(2)*y0)/sigma ) * tan ( (sigma*h)/sqrt(2) ) )
	m3 = s^3*exp(((sqrt(6)*y0)/sigma ) * tan ( (sigma*h*sqrt(6))/2 ) )
	m4 = s^4*exp(((sqrt(12)*y0)/sigma ) * tan (sigma*h*sqrt(3)) )

	c = (m1*m4-m2*m3)/(m1*m3-m2^2)
	q = (m2*m4-m3^2)/(m1*m3-m2^2)
	delta = c^2 - 4*q

	x1 = (c - sqrt(abs(delta)) ) / 2
	x2 = (c + sqrt(abs(delta)) ) / 2
	x3 = 0

	p1 = (m1*x2 - m2) / (x1*(x2-x1))

	p2 = (m2 - m1*x1) / (x2*(x2-x1))

	if (uniform < p1) {
		answer <- x1
	}
	else 
	if (uniform < p1+p2) {
		answer <- x2
	} else {
		answer <- x3
	}
	
	return (answer)
}

ad3 <- function (s, y0, h, sigma, uniform) {
# Senas bandymas 4 atsitiktiniais dydziais generuoti S_h
	#Naudojam 4 vidurkius...

	m1 = s
	m2 = s*s*exp(((sqrt(2)*y0)/sigma ) * tan ( (sigma*h)/sqrt(2) ) )
	m3 = s^3*exp(((sqrt(6)*y0)/sigma ) * tan ( (sigma*h*sqrt(6))/2 ) )
	m4 = s^4*exp(((sqrt(12)*y0)/sigma ) * tan (sigma*h*sqrt(3)) )

	y1 = 0
	y2 = m2 + s^2
	y3 = m3 - s^3-3*s*m2
	y4 = m4 + s^4 + 6*s^2*m2 - 4*s*m3

	c = (m3-s*(s^2+3*m2))/(m2+s^2)
	delta = (-1)*3*c^2 + 4*(y4)/(y2)

	x1_plus = (c - sqrt(abs(delta)) ) / 2
	x2_plus = (c + sqrt(abs(delta)) ) / 2
	x3_plus = 0

	p1 = ( (-1)*y2) / (x1_plus*(x2_plus-x1_plus))
	p2 = (y2) / (x2_plus*(x2_plus-x1_plus))

	print ("ad3")
	print (s)
	print (delta)
	print (x1_plus)
	print (x2_plus)
	print (p1)
	print (p2)

	if (uniform < p1) {
		answer <- (s+x1_plus)
	}
	else 
	if (uniform < p1+p2) {
		answer <- (s+x2_plus)
	} else {
		answer <- s
	}
	
	return (answer)
}

CIRD <- function (x, h, teta, k) {
  # Deterministines CIR dalies sprendimas
  
  e = exp ((-1) * h * k)
  # Toks atsakymas butu, jei CIR butu (teta - k X_t) dt:
  #answer <- (x * e + teta/k * (1-e) )
  # Musu atveju yra k(teta - X_t) dt:
  answer <- (x * e + teta * (1-e) )
  
  return (answer)
}

XD <- function (x, y, h, k, r, teta) {
  # Hestono modelio pirmosios lygties X_t = log(S_t) deterministines dalies sprendimas  
  
  e = exp ((-1) * h * k)
  
  answer <- (x + h * (r-teta/2) + (y - teta) / (2*k) * (e-1) )	
  
  #	print (answer)
  
  #if (answer < 0)
  #	print ("Neigiamas D")
  
  return (answer)
}



D <- function (x, h, r) {
# Hestono modelio pirmosios lygties (S_t) deterministines dalies sprendimas  

	answer <- (x * exp (h * r) )	

#	print (answer)

	#if (answer < 0)
	#	print ("Neigiamas D")

	return (answer)
}

tikra_funkcija <- function (s0, t, sigma) {
# OBSOLETE
	answer <- s0*s0* exp ( ( (sqrt(2)*s0)/sigma ) * tan ( (sigma*t)/sqrt(2) ) )
	return (answer)
}
tikras_sprendinys <- function (s0, x, sigma) {
# OBSOLETE	
	print (length(x))
	for (i in 1:(length(x)-1)) {
		y = c(y, tikra_funkcija (s0, x[i], sigma))
	}	
	
	return (y)
}

tikras_momentas2 <- function (s0, y0, t) {
# OBSOLETE
	psi = tan (atan(1)-t/2) - 1
	phi = (-2)*log(cos(atan(1)-t/2))-t

	answer = s0*s0*exp(t)*exp((-1)*phi)*exp((-1)*y0*psi)
	
	return (answer)
}

tikras_momentas_n_rho <- function (s0, y0, n, r, k, teta, sigma, rho, t) {
# Tikra n-tojo momento reiksme taske t.	
# rho <> 0
	
	chi = rho * sigma * n - k
	delta = chi^2 - sigma^2*(n^2-n)

	print ("___chi__delta_______")
	print (chi)
	print (delta)  
  
	if (delta < 0) {

		C0 = atan (chi/sqrt ( (-1) * delta))
		t_sigma = sqrt ( (-1) * delta) * t / 2
		phi = ( (2*k*teta)/(sigma^2) ) * log (abs (cos (C0))) - ((2*k*teta)/ (sigma^2)) * log (abs (cos(t_sigma + C0))) + (r*n-(chi*k*teta)/(sigma^2))*t
		psi = sqrt ( (-1) * delta)/(sigma^2) * tan (t_sigma + C0) - chi/(sigma^2)

	} else
	if (delta > 0) {
    
    C = (chi - sqrt(delta))/(chi + sqrt(delta))
    
		gama = log( C )
		koef = (k*teta)
    #e = exp (t*sqrt(delta))
    
    minuse = exp ((-1)*t*sqrt(delta))
    
    phi1 = (sqrt(delta) - chi) / (sigma^2)
    phi2 = (-1) * (sqrt(delta) + chi) / (sigma^2)
    #e = ((sigma^2) * (phi1-phi2) *t) / 2
    e = exp (t*sqrt(delta))
    
    
    #e = exp (t*sqrt(delta)+ gama)
		#phi = t * (r+1-koef) - 2*koef * log ((2*sqrt(delta) * (1 - exp(t*sqrt(delta)+ C0)))/(chi + sqrt(delta)))
		
    #phi = t * (r*n+ (koef*sqrt(delta)-koef*chi)/(sigma^2) ) - 2*(koef/sqrt(delta)) * log ( abs( (1 - e)*(chi + sqrt(delta)) ) / (2*sqrt(delta)) )
		
    #phi = t * (r*n+koef-(koef*chi/(sigma^2)) ) - 2*(koef / sqrt(delta)) * log ( abs( (1 + e)*(chi + sqrt(delta)) ) / 2*chi )

    #psi = (sqrt(delta) * (exp(t*sqrt(delta) + C0) + 1) ) / (sigma^2*(1-exp(t*sqrt(delta)+C0))) - chi/(sigma^2)
    
    #psi = (sqrt(delta) * (1+e)) / (sigma^2* (1-e)) - chi/(sigma^2)  
		
    #psi = (sqrt(delta) * (1-e)) / (sigma^2* (1+e)) - chi/(sigma^2)
    
    phi = t*(r*n - (koef*(chi+sqrt(delta))/sigma^2) ) - (2*koef)/(sigma^2) * log ( (sqrt(delta) - chi + (chi+sqrt(delta))*minuse )/( (2)*sqrt(delta) )  ) 
    #psi = ( (delta - chi^2)*(1-e) ) / ( sigma^2 * ( chi + sqrt(delta) - (chi - sqrt(delta))* e) )
    
    psi = phi2 + (phi2-phi1) / ((phi1/phi2) * e - 1)
    
    test = sigma^2 * psi + chi
    sqrtDelta = sqrt(delta)
    
    print ("----test-psi--chi--sqrt(Delta)------")
    print (test)
		print (psi)
    print (chi)
    print (sqrtDelta)
    print ("----------------------------")
        
	} else {
		phi = (-1) * (k*teta*chi^2*t) / (sigma^2*(t*chi+2)) + r*t
		psi = (-1) * (chi^2*t) / (sigma^2*(t*chi+2))
	}
	
	answer = s0^n * exp(phi) * exp(y0*psi)

	return (answer)
}

tikras_momentas_n <- function (s0, y0, n, r, k, teta, sigma, t) {
# Tikra n-tojo momento reiksme taske t.	

	mu = (1/2) * n * (n-1)
	alpha = (2*mu*sigma*sigma - k*k)/(sigma^4)
	
  print ("___MU__ALFA_______")
	print (mu)
  print (alpha)
  
	if (alpha > 0) {
		C0 = atan (k/(sigma^2*sqrt(alpha)))
		phi = 1/(sigma*sigma) * (2*log(cos(C0 - t*((sigma^2)*sqrt(alpha))/2)) - 2*log(cos(C0)) - k*t)
		psi = sqrt(alpha) * tan(C0 - t*((sigma^2)*sqrt(alpha))/2) - k/(sigma^2)

	} else 
	if (alpha < 0) {
		C0 = log (abs ( (k-sigma^2*sqrt(-alpha))/(k+sigma^2*sqrt(-alpha)) ) )
		phi = t*(sqrt(-alpha) - k/sigma^2) + 4/(sigma^2*sqrt(-alpha)) * log (abs ( (1 - exp(C0-1/2*sigma^2*sqrt(-alpha)*t) )/ (1-exp(C0)) ) )
		psi = sqrt(-alpha) * ( exp(C0-1/2*sigma^2*sqrt(-alpha)*t) + 1 )/( 1 - exp(C0-1/2*sigma^2*sqrt(-alpha)*t) ) - k/sigma^2
	}
	else {
		phi = 1/sigma^2 * (2* log ( abs ( (2+k*t)/2 )) - k*t)
		psi = (-k^2*t)/(sigma^2*(2+k*t))
	}

	answer = s0^n * exp(n*r*t) * exp((-1)*k*teta*phi) * exp((-1) * y0 * psi)
	return (answer)
}

tikras_momentas_log <- function (s0, y0, k, teta,  t) {
# Tikra vidurkio E ln (S_t) reiksme
	answer = log(s0) + 1/(2*k) * (teta - y0) * (1 - exp ((-1)*k*t) ) - (1/2) * teta * t 
	
	return (answer)
}






